In this notebook, we conduct additional analyses for Study 1: In-depth interviews with people of faith (not reported in the main text or the Supplemental Information of the paper).
source("../../scripts_general/dependencies.R")
source("../../scripts_general/custom_funs.R")
source("../../scripts_general/var_recode_contrast.R")
source("../scripts_s1/s1_var_groups.R")
source("../../scripts_general/data_load.R")
Parsed with column specification:
cols(
.default = col_double(),
date = [34mcol_date(format = "")[39m,
researcher = [31mcol_character()[39m,
country = [31mcol_character()[39m,
site = [31mcol_character()[39m,
religion = [31mcol_character()[39m,
subject_gender = [31mcol_character()[39m,
subject_job = [31mcol_character()[39m,
subject_schedule = [31mcol_character()[39m,
subject_livedhere = [31mcol_character()[39m,
subject_lang = [31mcol_character()[39m,
subject_marital = [31mcol_character()[39m,
subject_hs = [31mcol_character()[39m,
subject_liveswith = [31mcol_character()[39m,
servicesperweek = [31mcol_character()[39m,
study = [31mcol_character()[39m
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
.default = col_double(),
date = [34mcol_date(format = "")[39m,
researcher = [31mcol_character()[39m,
country = [31mcol_character()[39m,
quad = [31mcol_character()[39m,
subject_gender = [31mcol_character()[39m,
subject_job = [31mcol_character()[39m,
subject_schedule = [31mcol_character()[39m,
subject_livedhere = [31mcol_character()[39m,
subject_lang = [31mcol_character()[39m,
subject_marital = [31mcol_character()[39m,
subject_hs = [31mcol_character()[39m,
subject_school = [31mcol_character()[39m,
subject_liveswith = [31mcol_character()[39m,
servicesperweek = [31mcol_character()[39m,
godexpviaawe_freq = [31mcol_character()[39m,
sleephabit = [31mcol_character()[39m,
researcher_date = [34mcol_date(format = "")[39m,
researcher_notes = [31mcol_character()[39m,
site = [31mcol_character()[39m,
religion = [31mcol_character()[39m
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
.default = col_double(),
ctry = [31mcol_character()[39m,
wher = [31mcol_character()[39m,
recr = [31mcol_character()[39m,
whoc = [31mcol_character()[39m,
demo_chur = [31mcol_character()[39m,
demo_ethn = [31mcol_character()[39m,
demo_maj = [31mcol_character()[39m,
demo_pocc = [31mcol_character()[39m,
demo_rlgn = [31mcol_character()[39m,
demo_sex = [31mcol_character()[39m
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
.default = col_character(),
X1 = [32mcol_double()[39m,
epi_subj = [32mcol_double()[39m,
epi_demo_age = [32mcol_double()[39m,
epi_demo_ses_num = [32mcol_double()[39m,
epi_demo_howr_num = [32mcol_double()[39m,
epi_demo_affr_num = [32mcol_double()[39m,
epi_demo_tung_num = [32mcol_double()[39m,
epi_demo_ytng_num = [32mcol_double()[39m,
epi_demo_affr_cat = [33mcol_logical()[39m,
epi_demo_tung_cat = [33mcol_logical()[39m,
epi_demo_ytng_cat = [33mcol_logical()[39m,
epi_version = [32mcol_double()[39m,
epi_charc = [32mcol_double()[39m
)
See spec(...) for full column specifications.
2 parsing failures.
row col expected actual file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
X1 = [32mcol_double()[39m,
epi_ctry = [31mcol_character()[39m,
epi_subj = [32mcol_double()[39m,
score = [32mcol_double()[39m
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
X1 = [32mcol_double()[39m,
epi_ctry = [31mcol_character()[39m,
epi_subj = [32mcol_double()[39m,
score = [32mcol_double()[39m
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
X1 = [32mcol_double()[39m,
epi_ctry = [31mcol_character()[39m,
epi_subj = [32mcol_double()[39m,
question = [31mcol_character()[39m,
response = [32mcol_double()[39m,
order = [32mcol_double()[39m,
question_text = [31mcol_character()[39m
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
X1 = [32mcol_double()[39m,
epi_ctry = [31mcol_character()[39m,
epi_subj = [32mcol_double()[39m,
question = [31mcol_character()[39m,
response = [32mcol_double()[39m,
order = [32mcol_double()[39m,
question_text = [31mcol_character()[39m
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
.default = col_double(),
p7_ctry = [31mcol_character()[39m,
p7_abs_check = [31mcol_character()[39m,
p7_dse_check = [31mcol_character()[39m,
p7_se_check = [31mcol_character()[39m,
p7_unev_check = [31mcol_character()[39m,
p7_exsen_check = [31mcol_character()[39m,
p7_por_check = [31mcol_character()[39m,
p7_mm_check = [31mcol_character()[39m,
p7_dem_sex = [31mcol_character()[39m,
p7_dem_pocc = [31mcol_character()[39m,
p7_dem_major = [31mcol_character()[39m,
p7_dem_ethnicity = [31mcol_character()[39m,
p7_dem_rur.urb = [31mcol_character()[39m,
p7_dem_affrd.basics = [31mcol_character()[39m,
p7_dem_religion = [31mcol_character()[39m,
p7_dem_church = [31mcol_character()[39m,
p7_dem_holy.tung.gif = [31mcol_character()[39m,
p7_abs_child.exp_cat = [33mcol_logical()[39m,
p7_abs_poetic_cat = [33mcol_logical()[39m,
p7_abs_tv.real_cat = [33mcol_logical()[39m
# ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
d1r_byquestion <- read_csv("../data/study1r_byquestion.csv")
Parsed with column specification:
cols(
.default = col_double(),
researcher = [31mcol_character()[39m,
country = [31mcol_character()[39m,
quad = [31mcol_character()[39m,
religion = [31mcol_character()[39m,
coder = [31mcol_character()[39m,
site = [31mcol_character()[39m
)
See spec(...) for full column specifications.
d1r <- read_csv("../data/study1r.csv") %>%
rename_at(vars(contains("score")), funs(gsub("score", "recoded", .))) %>%
full_join(d1) %>%
mutate(country = factor(country, levels = levels_country),
site = factor(site, levels = levels_site),
religion = factor(religion, levels = levels_religion),
researcher = factor(researcher, levels = levels_researcher))
Parsed with column specification:
cols(
researcher = [31mcol_character()[39m,
country = [31mcol_character()[39m,
site = [31mcol_character()[39m,
religion = [31mcol_character()[39m,
coder = [31mcol_character()[39m,
subject_id = [32mcol_double()[39m,
spirit_score = [32mcol_double()[39m,
other_score = [32mcol_double()[39m,
study = [31mcol_character()[39m,
spirit_score_std = [32mcol_double()[39m,
other_score_std = [32mcol_double()[39m,
spirit_score_std2 = [32mcol_double()[39m,
other_score_std2 = [32mcol_double()[39m
)
Joining, by = c("researcher", "country", "site", "religion", "subject_id", "study")
Column `researcher` joining character vector and factor, coercing into character vectorColumn `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector
contrasts(d1r$country) <- contrasts_country
contrasts(d1r$site) <- contrasts_site
contrasts(d1r$religion) <- contrasts_religion
d1r %>%
select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
gather(scale, score, -c(country, subject_id)) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
grepl("other", scale) ~ "other extraordinary experiences",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual experiences",
"other extraordinary experiences"))) %>%
ggplot(aes(x = scale, y = score, group = score_type)) +
# geom_point(aes(color = score_type), alpha = 0.2,
# position = position_jitterdodge(jitter.width = 0.1,
# dodge.width = 0.5,
# jitter.height = 0.04)) +
geom_pointrange(data = . %>%
group_by(scale, score_type) %>%
langcog::multi_boot_standard(col = "score", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
position = position_dodge(width = 0.5),
shape = 21, fatten = 6) +
labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
x = "Scale", y = "Score (range: 0-1)",
fill = "Score type")
lmer(score ~ score_type * scale + (1 | subject_id),
data = d1r %>%
select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
gather(scale, score, -c(country, subject_id)) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
grepl("other", scale) ~ "other extraordinary experiences",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual experiences",
"other extraordinary experiences"))),
contrasts = list(score_type = "contr.sum", scale = "contr.sum")) %>%
summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: score ~ score_type * scale + (1 | subject_id)
Data:
d1r %>% select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
gather(scale, score, -c(country, subject_id)) %>% filter(grepl("spirit",
scale) | grepl("other", scale)) %>% mutate(score_type = case_when(grepl("recoded",
scale) ~ "recoded", grepl("score", scale) ~ "original", TRUE ~
NA_character_), score_type = factor(score_type, levels = c("original",
"recoded")), scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
grepl("other", scale) ~ "other extraordinary experiences",
TRUE ~ NA_character_), scale = factor(scale, levels = c("spiritual experiences",
"other extraordinary experiences")))
REML criterion at convergence: 362.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.37869 -0.68474 0.01642 0.55450 2.82615
Random effects:
Groups Name Variance Std.Dev.
subject_id (Intercept) 0.04320 0.2078
Residual 0.05204 0.2281
Number of obs: 1336, groups: subject_id, 341
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.251e-01 1.289e-02 3.416e+02 32.992 < 2e-16 ***
score_type1 -4.429e-03 6.269e-03 1.003e+03 -0.706 0.480112
scale1 2.433e-02 6.252e-03 9.972e+02 3.892 0.000106 ***
score_type1:scale1 6.545e-03 6.246e-03 9.954e+02 1.048 0.294921
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) scr_t1 scale1
score_type1 0.004
scale1 -0.006 0.006
scr_typ1:s1 0.003 -0.010 0.003
d1r %>%
select(country, subject_id, ends_with("recoded"), ends_with("score")) %>%
gather(scale, score, -c(country, subject_id)) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
grepl("other", scale) ~ "other exp.",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual exp.",
"other exp."))) %>%
ggplot(aes(x = scale, y = score, group = score_type)) +
facet_grid(~ country) +
# geom_point(aes(color = score_type), alpha = 0.2,
# position = position_jitterdodge(jitter.width = 0.1,
# dodge.width = 0.5,
# jitter.height = 0.04)) +
geom_pointrange(data = . %>%
group_by(scale, score_type, country) %>%
langcog::multi_boot_standard(col = "score", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
position = position_dodge(width = 0.5),
shape = 21, fatten = 6) +
labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
x = "Scale", y = "Score (range: 0-1)",
fill = "Score type") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
d1r %>%
select(country, site, religion, subject_id,
ends_with("recoded"), ends_with("score")) %>%
gather(scale, score, -c(country, site, religion, subject_id)) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
grepl("other", scale) ~ "other exp.",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual exp.",
"other exp."))) %>%
ggplot(aes(x = scale, y = score, group = score_type)) +
facet_grid(rows = vars(site, religion), cols = vars(country)) +
geom_point(aes(color = score_type), alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
dodge.width = 0.5,
jitter.height = 0.04)) +
geom_pointrange(data = . %>%
group_by(scale, score_type, country, site, religion) %>%
langcog::multi_boot_standard(col = "score", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = score_type),
position = position_dodge(width = 0.5),
shape = 21, fatten = 6) +
labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
x = "Scale", y = "Score (range: 0-1)",
fill = "Score type", color = "Score type") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
d1r %>%
select(country, subject_id, coder, ends_with("recoded")) %>%
gather(scale, score, -c(country, subject_id, coder)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
grepl("other", scale) ~ "other exp.",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual exp.",
"other exp."))) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
ggplot(aes(x = scale, y = score, group = coder)) +
facet_grid(~ country) +
geom_point(aes(color = coder), alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
dodge.width = 0.5,
jitter.height = 0.04)) +
geom_pointrange(data = . %>%
group_by(scale, coder, country) %>%
langcog::multi_boot_standard(col = "score", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = coder),
shape = 21, fatten = 5,
position = position_dodge(width = 0.5)) +
labs(title = "Comparing coders in recoded dataset (QTS2)",
x = "Scale", y = "Score (range: 0-1)", fill = "Coder", color = "Coder") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
lmer(score ~ scale * coder + (1 | subject_id),
data = d1r %>%
select(country, coder, subject_id, ends_with("recoded")) %>%
distinct() %>%
gather(scale, score, -c(country, subject_id, coder)) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
mutate(scale = case_when(grepl("spirit", scale) ~ "spiritual experiences",
grepl("other", scale) ~ "other extraordinary experiences",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual experiences",
"other extraordinary experiences")),
coder = factor(coder, levels = c("Nikki", "Maria", "Lucy"))),
contrasts = list(scale = "contr.sum",
coder = cbind("_NvGM" = c(1, 0, -1),
"_MvGM" = c(0, 1, -1)))) %>%
summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: score ~ scale * coder + (1 | subject_id)
Data: d1r %>% select(country, coder, subject_id, ends_with("recoded")) %>%
distinct() %>% gather(scale, score, -c(country, subject_id,
coder)) %>% filter(grepl("spirit", scale) | grepl("other",
scale)) %>% mutate(scale = case_when(grepl("spirit", scale) ~
"spiritual experiences", grepl("other", scale) ~ "other extraordinary experiences",
TRUE ~ NA_character_), scale = factor(scale, levels = c("spiritual experiences",
"other extraordinary experiences")), coder = factor(coder,
levels = c("Nikki", "Maria", "Lucy")))
REML criterion at convergence: 328.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.76292 -0.68127 0.03902 0.47704 2.27078
Random effects:
Groups Name Variance Std.Dev.
subject_id (Intercept) 0.02397 0.1548
Residual 0.06998 0.2645
Number of obs: 671, groups: subject_id, 340
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.42875 0.01515 341.70286 28.298 < 2e-16 ***
scale1 0.02745 0.01174 339.27621 2.338 0.01995 *
coder_NvGM 0.02927 0.01839 340.24072 1.591 0.11246
coder_MvGM -0.06042 0.02038 340.51182 -2.965 0.00324 **
scale1:coder_NvGM -0.03378 0.01423 337.85385 -2.374 0.01818 *
scale1:coder_MvGM 0.01857 0.01577 338.11808 1.178 0.23978
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) scale1 cd_NGM cd_MGM s1:_NG
scale1 -0.022
coder_NvGM -0.434 0.015
coder_MvGM -0.142 0.008 -0.173
scl1:cd_NGM 0.015 -0.437 -0.017 -0.004
scl1:cd_MGM 0.008 -0.145 -0.004 -0.018 -0.169
d1r %>%
select(country, site, religion, subject_id, coder, ends_with("recoded")) %>%
gather(scale, score, -c(country, site, religion, subject_id, coder)) %>%
mutate(score_type = case_when(grepl("recoded", scale) ~ "recoded",
grepl("score", scale) ~ "original",
TRUE ~ NA_character_),
score_type = factor(score_type, levels = c("original", "recoded")),
scale = case_when(grepl("spirit", scale) ~ "spiritual exp.",
grepl("other", scale) ~ "other exp.",
TRUE ~ NA_character_),
scale = factor(scale, levels = c("spiritual exp.",
"other exp."))) %>%
filter(grepl("spirit", scale) | grepl("other", scale)) %>%
ggplot(aes(x = scale, y = score, group = coder)) +
facet_grid(rows = vars(site, religion), cols = vars(country)) +
geom_point(aes(color = coder), alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
dodge.width = 0.5,
jitter.height = 0.04)) +
geom_pointrange(data = . %>%
group_by(scale, coder, country, site, religion) %>%
langcog::multi_boot_standard(col = "score", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper, fill = coder),
shape = 21, fatten = 5,
position = position_dodge(width = 0.5)) +
geom_text(data = . %>% count(country, site, religion, coder, scale),
aes(y = 0, label = paste0("n=", n), color = coder),
position = position_dodge(width = 1)) +
labs(title = "Comparing coders in recoded dataset (QTS2)",
x = "Scale", y = "Score (range: 0-1)", fill = "Coder", color = "Coder") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
d1r_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, recoded, -c(country, site, religion, subject_id)) %>%
full_join(d1_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, original, -c(country, site, religion, subject_id))) %>%
gather(response_type, response, c(original, recoded)) %>%
mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other)),
response_type = factor(response_type, levels = c("original", "recoded"))) %>%
ggplot(aes(x = question, y = response, group = response_type)) +
geom_pointrange(data = . %>%
group_by(question, response_type) %>%
langcog::multi_boot_standard(col = "response", na.rm = T),
aes(y = mean, ymin = ci_lower, ymax = ci_upper,
fill = response_type),
position = position_dodge(width = 0.5),
shape = 21, fatten = 6) +
labs(title = "Comparing recoded scores (QTS2) to original scores (QTS1)",
x = "Question", y = "Response (range: 0-1)",
fill = "Response type") +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector
d1r_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, recoded, -c(country, site, religion, subject_id)) %>%
full_join(d1_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, original, -c(country, site, religion, subject_id))) %>%
gather(response_type, response, c(original, recoded)) %>%
mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other)),
response_type = factor(response_type, levels = c("original", "recoded"))) %>%
filter(question == "seehearnotgod") %>%
count(response_type, response) %>%
spread(response, n)
Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector
d1r_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, recoded, -c(country, site, religion, subject_id)) %>%
full_join(d1_byquestion %>%
select(country, site, religion, subject_id,
!!s1_var_spirit, !!s1_var_other) %>%
gather(question, original, -c(country, site, religion, subject_id))) %>%
mutate(question = factor(question, levels = c(s1_var_spirit, s1_var_other))) %>%
filter(question == "seehearnotgod") %>%
mutate(match = (original == recoded)) %>%
count(match) %>%
mutate(prop = n/sum(n))
Joining, by = c("country", "site", "religion", "subject_id", "question")
Column `country` joining character vector and factor, coercing into character vectorColumn `site` joining character vector and factor, coercing into character vectorColumn `religion` joining character vector and factor, coercing into character vector
d1_check_questions <- c("godvoxaloud", "godviavisions", "godviatouch",
"godviasmell", "neartangiblegod", "presencenotgod",
"presencedemon", "beingentbody", "seehearnotgod",
"whitelight", "spiritbeingencounter", "voxwhenalone",
"seethingscornereye")
d1_check <- d1_byquestion %>%
select(country, subject_id, ends_with("_judge")) %>%
gather(question, judgment, -c(country, subject_id)) %>%
mutate(question = gsub("_judge", "", question)) %>%
full_join(d1_byquestion %>%
select(country, subject_id, !!d1_check_questions) %>%
gather(question, response, -c(country, subject_id)))
Joining, by = c("country", "subject_id", "question")
d1_check %>%
filter(question %in% d1_check_questions) %>%
mutate(agree = (judgment == response),
agree = ifelse(is.na(agree), 0, agree)) %>%
group_by(question) %>%
summarise(percent_agree = mean(agree, na.rm = T)) %>%
arrange(percent_agree)
s1_check_icc_df <- data.frame(question = character(), icc = numeric())
for(i in d1_check_questions){
temp_icc <- icc_fun(d1_check, var_name = i,
var1 = "judgment", var2 = "response")
temp_res <- data.frame(question = i, icc = temp_icc)
s1_check_icc_df <- full_join(s1_check_icc_df, temp_res)
rm(temp_res, i)
}
Joining, by = c("question", "icc")
Column `question` joining factors with different levels, coercing to character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vectorJoining, by = c("question", "icc")
Column `question` joining character vector and factor, coercing into character vector
s1_check_icc_df <- s1_check_icc_df %>%
mutate(question = factor(question, levels = d1_check_questions),
reliability = case_when(icc < 0.5 ~ "poor",
icc < 0.75 ~ "moderate",
icc < 0.9 ~ "good",
icc <= 1 ~ "excellent",
TRUE ~ NA_character_),
reliability = factor(reliability,
levels = c("poor", "moderate",
"good", "excellent")),
scale = case_when(question %in% s1_var_spirit ~ "spiritual experiences",
question %in% s1_var_other ~ "other experiences",
TRUE ~ NA_character_),
scale = factor(scale,
levels = c("spiritual experiences", "other experiences")))
s1_check_icc_df %>% arrange(icc)
s1_check_icc_df %>%
mutate(scale = recode_factor(scale,
"spiritual experiences" = "spiritual",
"other experiences" = "other")) %>%
ggplot(aes(x = question, y = icc, color = reliability)) +
facet_grid(~ scale, scales = "free", space = "free") +
geom_point() +
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2)) +
scale_color_brewer(palette = "RdYlGn", direction = 1, drop = F) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1,
color = ifelse(s1_check_icc_df$reliability == "poor",
"#d7191c", "black"))) +
labs(title = "Comparing initial responses vs. later 'judgments'",
x = "Question", y = "Intraclass correlation coefficient (ICC)")
d1_byquestion %>%
count(godvoxaloud) %>%
mutate(godvoxaloud = recode_factor(godvoxaloud,
"0" = "no", "0.5" = "maybe", "1" = "yes"),
proportion = round(n/sum(n), 2))
d1_byquestion %>%
filter(!is.na(godvoxaloud)) %>%
count(country, godvoxaloud) %>%
mutate(godvoxaloud = recode_factor(godvoxaloud,
"0" = "no", "0.5" = "maybe", "1" = "yes")) %>%
complete(country, nesting(godvoxaloud), fill = list(n = 0)) %>%
spread(godvoxaloud, n) %>%
mutate(total = no + maybe + yes) %>%
janitor::adorn_totals() %>%
mutate_at(vars(no, maybe, yes), funs(round(./total, 2)*100)) %>%
select(-total) %>%
data.frame() %>%
mutate_at(vars(no, maybe, yes), funs(paste0(., "%"))) %>%
kable() %>%
kable_styling() %>%
row_spec(6, bold = T)
Column `country` has different attributes on LHS and RHS of join
| country | no | maybe | yes |
|---|---|---|---|
| US | 86% | 1% | 13% |
| Ghana | 33% | 0% | 67% |
| Thailand | 58% | 2% | 40% |
| China | 50% | 3% | 47% |
| Vanuatu | 54% | 9% | 37% |
| Total | 57% | 3% | 41% |
d1_byquestion %>%
filter(!is.na(godviavisions)) %>%
count(country, godviavisions) %>%
mutate(godviavisions = recode_factor(godviavisions,
"0" = "no", "0.5" = "maybe", "1" = "yes")) %>%
complete(country, nesting(godviavisions), fill = list(n = 0)) %>%
spread(godviavisions, n) %>%
mutate(total = no + maybe + yes) %>%
janitor::adorn_totals() %>%
mutate_at(vars(no, maybe, yes), funs(round(./total, 2)*100)) %>%
select(-total) %>%
data.frame() %>%
mutate_at(vars(no, maybe, yes), funs(paste0(., "%"))) %>%
kable() %>%
kable_styling() %>%
row_spec(6, bold = T)
Column `country` has different attributes on LHS and RHS of join
| country | no | maybe | yes |
|---|---|---|---|
| US | 78% | 0% | 22% |
| Ghana | 41% | 0% | 59% |
| Thailand | 55% | 2% | 43% |
| China | 78% | 0% | 22% |
| Vanuatu | 61% | 0% | 39% |
| Total | 63% | 0% | 37% |
godvoxaloud and godviavisionsd1_tml1 <- d1 %>% select(country, subject_id, abs_score, por_score) %>%
full_join(d1_byquestion %>%
select(country, subject_id, godvoxaloud, godviavisions) %>%
gather(question, response, -c(country, subject_id)) %>%
mutate(response_yn = case_when(response < 1 ~ 0,
response == 1 ~ 1,
TRUE ~ NA_real_)))
Joining, by = c("country", "subject_id")
glmer(response_yn ~ scale(por_score) + (1 | question) + (1 | country),
family = "binomial", d1_tml1) %>% summary()
boundary (singular) fit: see ?isSingular
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(por_score) + (1 | question) + (1 | country)
Data: d1_tml1
AIC BIC logLik deviance df.resid
769.4 787.2 -380.7 761.4 620
Scaled residuals:
Min 1Q Median 3Q Max
-1.7952 -0.7330 -0.4624 0.8706 2.4182
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.297 0.545
question (Intercept) 0.000 0.000
Number of obs: 624, groups: country, 5; question, 2
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5236 0.2595 -2.018 0.04360 *
scale(por_score) 0.3109 0.1090 2.851 0.00436 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(pr_scr) -0.006
convergence code: 0
boundary (singular) fit: see ?isSingular
glmer(response_yn ~ scale(abs_score) + (1 | question) + (1 | country),
family = "binomial", d1_tml1) %>% summary()
boundary (singular) fit: see ?isSingular
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | question) + (1 | country)
Data: d1_tml1
AIC BIC logLik deviance df.resid
745.3 762.9 -368.6 737.3 603
Scaled residuals:
Min 1Q Median 3Q Max
-1.6946 -0.7569 -0.4474 0.9293 2.7417
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.4476 0.669
question (Intercept) 0.0000 0.000
Number of obs: 607, groups: country, 5; question, 2
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.53208 0.31285 -1.701 0.088989 .
scale(abs_score) 0.33016 0.09554 3.456 0.000549 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(bs_scr) -0.028
convergence code: 0
boundary (singular) fit: see ?isSingular
glmer(response_yn ~ scale(por_score) + (1 | country),
family = "binomial", d1_tml1 %>% filter(question == "godvoxaloud")) %>%
summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(por_score) + (1 | country)
Data: d1_tml1 %>% filter(question == "godvoxaloud")
AIC BIC logLik deviance df.resid
385.8 397.0 -189.9 379.8 310
Scaled residuals:
Min 1Q Median 3Q Max
-1.7575 -0.7684 -0.4053 0.7893 2.5614
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.585 0.7649
Number of obs: 313, groups: country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4694 0.3651 -1.286 0.199
scale(por_score) 0.2172 0.1551 1.401 0.161
Correlation of Fixed Effects:
(Intr)
scl(pr_scr) 0.000
glmer(response_yn ~ scale(abs_score) + (1 | country),
family = "binomial", d1_tml1 %>% filter(question == "godvoxaloud")) %>%
summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | country)
Data: d1_tml1 %>% filter(question == "godvoxaloud")
AIC BIC logLik deviance df.resid
375.5 386.7 -184.8 369.5 302
Scaled residuals:
Min 1Q Median 3Q Max
-1.6176 -0.7766 -0.3915 0.8313 2.8594
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.7087 0.8418
Number of obs: 305, groups: country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4886 0.3986 -1.226 0.220
scale(abs_score) 0.2052 0.1343 1.529 0.126
Correlation of Fixed Effects:
(Intr)
scl(bs_scr) -0.019
glmer(response_yn ~ scale(por_score) + (1 | country),
family = "binomial", d1_tml1 %>% filter(question == "godviavisions")) %>%
summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(por_score) + (1 | country)
Data: d1_tml1 %>% filter(question == "godviavisions")
AIC BIC logLik deviance df.resid
384.3 395.5 -189.1 378.3 308
Scaled residuals:
Min 1Q Median 3Q Max
-1.6899 -0.7062 -0.5004 0.9736 2.2955
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.09818 0.3133
Number of obs: 311, groups: country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6004 0.1891 -3.175 0.00150 **
scale(por_score) 0.5138 0.1657 3.100 0.00194 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(pr_scr) -0.017
glmer(response_yn ~ scale(abs_score) + (1 | country),
family = "binomial", d1_tml1 %>% filter(question == "godviavisions")) %>%
summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | country)
Data: d1_tml1 %>% filter(question == "godviavisions")
AIC BIC logLik deviance df.resid
371.7 382.9 -182.9 365.7 299
Scaled residuals:
Min 1Q Median 3Q Max
-1.6618 -0.7152 -0.4837 0.9376 2.9033
Random effects:
Groups Name Variance Std.Dev.
country (Intercept) 0.3431 0.5857
Number of obs: 302, groups: country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6102 0.2933 -2.081 0.037477 *
scale(abs_score) 0.4761 0.1387 3.432 0.000599 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(bs_scr) -0.066
d1_tml2 <- d1 %>% select(country, subject_id, abs_score, por_score) %>%
full_join(d1_byquestion %>%
select(country, subject_id, neartangiblegod, presencenotgod,
presencedemon, beingentbody, seehearnotgod,
spiritbeingencounter) %>%
gather(question, response, -c(country, subject_id)) %>%
mutate(response_yn = case_when(response < 1 ~ 0,
response == 1 ~ 1,
TRUE ~ NA_real_)))
Joining, by = c("country", "subject_id")
glmer(response_yn ~ scale(por_score) + (1 | question) + (1 | country),
family = "binomial", d1_tml2) %>% summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(por_score) + (1 | question) + (1 | country)
Data: d1_tml2
AIC BIC logLik deviance df.resid
2025.9 2048.1 -1009.0 2017.9 1875
Scaled residuals:
Min 1Q Median 3Q Max
-2.7229 -0.6165 -0.3640 0.7589 4.0868
Random effects:
Groups Name Variance Std.Dev.
question (Intercept) 0.4731 0.6878
country (Intercept) 0.6302 0.7939
Number of obs: 1879, groups: question, 6; country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.75940 0.45630 -1.664 0.0961 .
scale(por_score) 0.44752 0.06785 6.596 4.22e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(pr_scr) -0.012
glmer(response_yn ~ scale(abs_score) + (1 | question) + (1 | country),
family = "binomial", d1_tml2) %>% summary()
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
Family: binomial ( logit )
Formula: response_yn ~ scale(abs_score) + (1 | question) + (1 | country)
Data: d1_tml2
AIC BIC logLik deviance df.resid
1978.2 2000.3 -985.1 1970.2 1822
Scaled residuals:
Min 1Q Median 3Q Max
-2.6347 -0.6291 -0.3769 0.7575 4.4684
Random effects:
Groups Name Variance Std.Dev.
question (Intercept) 0.4829 0.6949
country (Intercept) 0.6679 0.8172
Number of obs: 1826, groups: question, 6; country, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.79078 0.46635 -1.696 0.0899 .
scale(abs_score) 0.28715 0.05964 4.814 1.48e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
scl(bs_scr) -0.015
d1_variance <- d1 %>%
select(country, site, religion, subject_id, abs_score, por_score) %>%
gather(var, score, ends_with("_score")) %>%
filter(!is.na(score)) %>%
mutate(score_rescaled = case_when(var == "abs_score" ~ score/1,
var == "por_score" ~ score/3)) # called pv_score in study 4
bartlett.test(score_rescaled ~ var,
data = d1_variance %>% filter(country == "US"))
Bartlett test of homogeneity of variances
data: score_rescaled by var
Bartlett's K-squared = 2.3384, df = 1, p-value = 0.1262
d1_variance %>%
filter(country == "US") %>%
group_by(var) %>%
summarise(variance = var(score_rescaled))
bartlett.test(score_rescaled ~ var,
data = d1_variance %>% filter(country == "Ghana"))
Bartlett test of homogeneity of variances
data: score_rescaled by var
Bartlett's K-squared = 4.6186, df = 1, p-value = 0.03163
d1_variance %>%
filter(country == "Ghana") %>%
group_by(var) %>%
summarise(variance = var(score_rescaled))
bartlett.test(score_rescaled ~ var,
data = d1_variance %>% filter(country == "Thailand"))
Bartlett test of homogeneity of variances
data: score_rescaled by var
Bartlett's K-squared = 0.8459, df = 1, p-value = 0.3577
d1_variance %>%
filter(country == "Thailand") %>%
group_by(var) %>%
summarise(variance = var(score_rescaled))
bartlett.test(score_rescaled ~ var,
data = d1_variance %>% filter(country == "China"))
Bartlett test of homogeneity of variances
data: score_rescaled by var
Bartlett's K-squared = 9.3242, df = 1, p-value = 0.002262
d1_variance %>%
filter(country == "China") %>%
group_by(var) %>%
summarise(variance = var(score_rescaled))
bartlett.test(score_rescaled ~ var,
data = d1_variance %>% filter(country == "Vanuatu"))
Bartlett test of homogeneity of variances
data: score_rescaled by var
Bartlett's K-squared = 16.04, df = 1, p-value = 6.202e-05
d1_variance %>%
filter(country == "Vanuatu") %>%
group_by(var) %>%
summarise(variance = var(score_rescaled))
d1_variance %>%
group_by(var, country) %>%
summarise(variance = var(score_rescaled)) %>%
ungroup() %>%
spread(var, variance) %>%
rename(Country = country, Absorption = abs_score,
`Porosity Vignettes` = por_score) %>%
kable(digits = 2) %>%
kable_styling()
| Country | Absorption | Porosity Vignettes |
|---|---|---|
| US | 0.04 | 0.03 |
| Ghana | 0.05 | 0.03 |
| Thailand | 0.03 | 0.02 |
| China | 0.05 | 0.02 |
| Vanuatu | 0.03 | 0.01 |
d1_variance %>%
group_by(var, country) %>%
mutate(score_cent = scale(score_rescaled, scale = F)) %>%
ungroup() %>%
mutate(var = recode_factor(var,
"por_score" = "Porosity Vignettes",
"abs_score" = "Absorption")) %>%
ggplot(aes(x = score_cent, fill = var, color = var, lty = var)) +
facet_wrap(country ~ ., ncol = 3) +
geom_density(alpha = 0.1) +
scale_linetype_manual(values = c(1, 2)) +
scale_color_brewer(palette = "Set1", direction = -1) +
scale_fill_brewer(palette = "Set1", direction = -1) +
labs(x = "Score (rescaled to run from 0-1, centered at the mean by country)",
color = "Scale", fill = "Scale", lty = "Scale") +
theme(legend.position = "bottom")
d1_variance %>%
group_by(var, country, site, religion) %>%
mutate(score_cent = scale(score_rescaled, scale = F)) %>%
ungroup() %>%
mutate(var = recode_factor(var,
"por_score" = "Porosity Vignettes",
"abs_score" = "Absorption")) %>%
ggplot(aes(x = score_cent, fill = var, color = var, lty = var)) +
facet_grid(cols = vars(country), rows = vars(site, religion)) +
geom_density(alpha = 0.1) +
scale_linetype_manual(values = c(1, 2)) +
scale_color_brewer(palette = "Set1", direction = -1) +
scale_fill_brewer(palette = "Set1", direction = -1) +
labs(x = "Score (rescaled to run from 0-1, centered at the mean by country)",
color = "Scale", fill = "Scale", lty = "Scale") +
theme(legend.position = "bottom")
r7 <- lm(por_score_std ~ country, d1)
regtab_fun(r7, country_var1 = "country_gh", country_var2 = "country_th", country_var3 = "country_ch", country_var4 = "country_vt")
rsquared(r7)
r9 <- lm(abs_score_std ~ country, d1)
rsquared(r9)
r10 <- lm(spirit_score_std ~ por_score_std * abs_score_std *
country * site * religion, d1)
# summary(r10)
regtab_fun(r10) %>% regtab_style_fun(row_emph = 2:3)
| Parameter | β | Std. Err. | t | p | |
|---|---|---|---|---|---|
| Intercept | 0.00 | 0.09 | -0.02 | 0.986 | |
| Porosity | 0.14 | 0.11 | 1.22 | 0.223 | |
| Absorption | 0.33 | 0.13 | 2.61 | 0.010 | ** |
| Country (Gh.) | -0.09 | 0.17 | -0.52 | 0.600 | |
| Country (Th.) | -0.03 | 0.20 | -0.17 | 0.867 | |
| Country (Ch.) | -0.39 | 0.23 | -1.68 | 0.095 | |
| Country (Vt.) | 0.86 | 0.18 | 4.76 | <0.001 | *** |
| Site (rural) | 0.09 | 0.09 | 0.95 | 0.345 | |
| Religion (Ch.) | 0.37 | 0.09 | 3.95 | <0.001 | *** |
| Porosity × Absorption | 0.00 | 0.19 | -0.02 | 0.987 | |
| Porosity × Country (Gh.) | 0.16 | 0.16 | 1.02 | 0.307 | |
| Porosity × Country (Th.) | 0.14 | 0.24 | 0.59 | 0.554 | |
| Porosity × Country (Ch.) | -0.17 | 0.22 | -0.76 | 0.446 | |
| Porosity × Country (Vt.) | -0.22 | 0.32 | -0.70 | 0.482 | |
| Absorption × Country (Gh.) | 0.17 | 0.19 | 0.91 | 0.365 | |
| Absorption × Country (Th.) | 0.23 | 0.29 | 0.79 | 0.429 | |
| Absorption × Country (Ch.) | 0.16 | 0.36 | 0.45 | 0.651 | |
| Absorption × Country (Vt.) | -0.33 | 0.21 | -1.63 | 0.106 | |
| Porosity × Site (rural) | -0.40 | 0.11 | -3.54 | <0.001 | *** |
| Absorption × Site (rural) | 0.00 | 0.13 | -0.03 | 0.978 | |
| Country (Gh.) × Site (rural) | -0.10 | 0.17 | -0.58 | 0.560 | |
| Country (Th.) × Site (rural) | -0.08 | 0.20 | -0.42 | 0.672 | |
| Country (Ch.) × Site (rural) | 0.22 | 0.23 | 0.94 | 0.350 | |
| Country (Vt.) × Site (rural) | 0.21 | 0.18 | 1.18 | 0.237 | |
| Porosity × Religion (Ch.) | -0.29 | 0.11 | -2.54 | 0.012 | * |
| Absorption × Religion (Ch.) | -0.23 | 0.13 | -1.79 | 0.075 | |
| Country (Gh.) × Religion (Ch.) | 0.04 | 0.17 | 0.22 | 0.824 | |
| Country (Th.) × Religion (Ch.) | -0.01 | 0.20 | -0.06 | 0.954 | |
| Country (Ch.) × Religion (Ch.) | -0.13 | 0.23 | -0.58 | 0.564 | |
| Country (Vt.) × Religion (Ch.) | 0.00 | 0.18 | 0.01 | 0.991 | |
| Site (rural) × Religion (Ch.) | -0.15 | 0.09 | -1.58 | 0.116 | |
| Porosity × Absorption × Country (Gh.) | -0.20 | 0.21 | -0.96 | 0.341 | |
| Porosity × Absorption × Country (Th.) | -0.51 | 0.51 | -1.00 | 0.316 | |
| Porosity × Absorption × Country (Ch.) | 0.18 | 0.43 | 0.42 | 0.673 | |
| Porosity × Absorption × Country (Vt.) | 0.53 | 0.41 | 1.31 | 0.192 | |
| Porosity × Absorption × Site (rural) | 0.20 | 0.19 | 1.09 | 0.277 | |
| Porosity × Country (Gh.) × Site (rural) | 0.16 | 0.16 | 1.00 | 0.318 | |
| Porosity × Country (Th.) × Site (rural) | 0.17 | 0.24 | 0.69 | 0.489 | |
| Porosity × Country (Ch.) × Site (rural) | 0.18 | 0.22 | 0.82 | 0.414 | |
| Porosity × Country (Vt.) × Site (rural) | -0.88 | 0.32 | -2.80 | 0.006 | ** |
| Absorption × Country (Gh.) × Site (rural) | -0.21 | 0.19 | -1.12 | 0.265 | |
| Absorption × Country (Th.) × Site (rural) | 0.30 | 0.29 | 1.04 | 0.302 | |
| Absorption × Country (Ch.) × Site (rural) | 0.00 | 0.36 | 0.00 | 1.000 | |
| Absorption × Country (Vt.) × Site (rural) | -0.24 | 0.21 | -1.15 | 0.252 | |
| Porosity × Absorption × Religion (Ch.) | -0.05 | 0.19 | -0.29 | 0.771 | |
| Porosity × Country (Gh.) × Religion (Ch.) | -0.03 | 0.16 | -0.21 | 0.831 | |
| Porosity × Country (Th.) × Religion (Ch.) | -0.25 | 0.24 | -1.05 | 0.293 | |
| Porosity × Country (Ch.) × Religion (Ch.) | 0.36 | 0.22 | 1.65 | 0.100 | |
| Porosity × Country (Vt.) × Religion (Ch.) | -0.33 | 0.32 | -1.04 | 0.301 | |
| Absorption × Country (Gh.) × Religion (Ch.) | -0.02 | 0.19 | -0.13 | 0.900 | |
| Absorption × Country (Th.) × Religion (Ch.) | 0.04 | 0.29 | 0.14 | 0.889 | |
| Absorption × Country (Ch.) × Religion (Ch.) | -0.34 | 0.36 | -0.94 | 0.349 | |
| Absorption × Country (Vt.) × Religion (Ch.) | 0.07 | 0.21 | 0.35 | 0.726 | |
| Porosity × Site (rural) × Religion (Ch.) | -0.10 | 0.11 | -0.86 | 0.391 | |
| Absorption × Site (rural) × Religion (Ch.) | 0.05 | 0.13 | 0.38 | 0.703 | |
| Country (Gh.) × Site (rural) × Religion (Ch.) | -0.52 | 0.17 | -3.06 | 0.002 | ** |
| Country (Th.) × Site (rural) × Religion (Ch.) | 0.53 | 0.20 | 2.68 | 0.008 | ** |
| Country (Ch.) × Site (rural) × Religion (Ch.) | -0.19 | 0.23 | -0.82 | 0.412 | |
| Country (Vt.) × Site (rural) × Religion (Ch.) | 0.08 | 0.18 | 0.44 | 0.659 | |
| Porosity × Absorption × Country (Gh.) × Site (rural) | 0.12 | 0.21 | 0.54 | 0.587 | |
| Porosity × Absorption × Country (Th.) × Site (rural) | -0.15 | 0.51 | -0.29 | 0.770 | |
| Porosity × Absorption × Country (Ch.) × Site (rural) | -0.32 | 0.43 | -0.75 | 0.457 | |
| Porosity × Absorption × Country (Vt.) × Site (rural) | 0.39 | 0.41 | 0.97 | 0.334 | |
| Porosity × Absorption × Country (Gh.) × Religion (Ch.) | 0.29 | 0.21 | 1.38 | 0.168 | |
| Porosity × Absorption × Country (Th.) × Religion (Ch.) | -0.53 | 0.51 | -1.03 | 0.302 | |
| Porosity × Absorption × Country (Ch.) × Religion (Ch.) | -0.38 | 0.43 | -0.89 | 0.375 | |
| Porosity × Absorption × Country (Vt.) × Religion (Ch.) | 0.43 | 0.41 | 1.06 | 0.290 | |
| Porosity × Absorption × Site (rural) × Religion (Ch.) | -0.14 | 0.19 | -0.75 | 0.454 | |
| Porosity × Country (Gh.) × Site (rural) × Religion (Ch.) | 0.45 | 0.16 | 2.84 | 0.005 | ** |
| Porosity × Country (Th.) × Site (rural) × Religion (Ch.) | 0.10 | 0.24 | 0.41 | 0.686 | |
| Porosity × Country (Ch.) × Site (rural) × Religion (Ch.) | -0.08 | 0.22 | -0.36 | 0.719 | |
| Porosity × Country (Vt.) × Site (rural) × Religion (Ch.) | -0.59 | 0.32 | -1.86 | 0.065 | |
| Absorption × Country (Gh.) × Site (rural) × Religion (Ch.) | 0.01 | 0.19 | 0.05 | 0.959 | |
| Absorption × Country (Th.) × Site (rural) × Religion (Ch.) | -0.02 | 0.29 | -0.07 | 0.946 | |
| Absorption × Country (Ch.) × Site (rural) × Religion (Ch.) | -0.10 | 0.36 | -0.27 | 0.784 | |
| Absorption × Country (Vt.) × Site (rural) × Religion (Ch.) | 0.02 | 0.21 | 0.08 | 0.937 | |
| Porosity × Absorption × Country (Gh.) × Site (rural) × Religion (Ch.) | 0.08 | 0.21 | 0.39 | 0.694 | |
| Porosity × Absorption × Country (Th.) × Site (rural) × Religion (Ch.) | -1.04 | 0.51 | -2.03 | 0.044 | * |
| Porosity × Absorption × Country (Ch.) × Site (rural) × Religion (Ch.) | 0.16 | 0.43 | 0.37 | 0.713 | |
| Porosity × Absorption × Country (Vt.) × Site (rural) × Religion (Ch.) | 0.44 | 0.41 | 1.08 | 0.279 |
r11 <- lm(spirit_score_std ~ por_score_std * abs_score_std * country * religion, d1)
# summary(r11)
regtab_fun(r11) %>% regtab_style_fun(row_emph = 2:3)
| Parameter | β | Std. Err. | t | p | |
|---|---|---|---|---|---|
| Intercept | -0.06 | 0.06 | -1.02 | 0.308 | |
| Porosity | 0.30 | 0.08 | 3.65 | <0.001 | *** |
| Absorption | 0.34 | 0.08 | 4.35 | <0.001 | *** |
| Country (Gh.) | 0.05 | 0.12 | 0.46 | 0.648 | |
| Country (Th.) | -0.04 | 0.11 | -0.32 | 0.750 | |
| Country (Ch.) | -0.39 | 0.16 | -2.45 | 0.015 | * |
| Country (Vt.) | 0.64 | 0.13 | 5.02 | <0.001 | *** |
| Religion (Ch.) | 0.31 | 0.06 | 4.93 | <0.001 | *** |
| Porosity × Absorption | 0.07 | 0.09 | 0.76 | 0.450 | |
| Porosity × Country (Gh.) | -0.14 | 0.12 | -1.21 | 0.226 | |
| Porosity × Country (Th.) | -0.01 | 0.15 | -0.07 | 0.944 | |
| Porosity × Country (Ch.) | -0.25 | 0.16 | -1.59 | 0.114 | |
| Porosity × Country (Vt.) | 0.44 | 0.24 | 1.84 | 0.067 | |
| Absorption × Country (Gh.) | -0.13 | 0.11 | -1.10 | 0.271 | |
| Absorption × Country (Th.) | 0.05 | 0.14 | 0.37 | 0.709 | |
| Absorption × Country (Ch.) | 0.39 | 0.23 | 1.67 | 0.095 | |
| Absorption × Country (Vt.) | -0.07 | 0.13 | -0.56 | 0.573 | |
| Porosity × Religion (Ch.) | -0.09 | 0.08 | -1.06 | 0.290 | |
| Absorption × Religion (Ch.) | -0.30 | 0.08 | -3.87 | <0.001 | *** |
| Country (Gh.) × Religion (Ch.) | -0.16 | 0.12 | -1.29 | 0.197 | |
| Country (Th.) × Religion (Ch.) | -0.22 | 0.11 | -1.99 | 0.047 | * |
| Country (Ch.) × Religion (Ch.) | 0.08 | 0.16 | 0.52 | 0.605 | |
| Country (Vt.) × Religion (Ch.) | 0.13 | 0.13 | 1.01 | 0.312 | |
| Porosity × Absorption × Country (Gh.) | -0.17 | 0.11 | -1.54 | 0.125 | |
| Porosity × Absorption × Country (Th.) | -0.27 | 0.18 | -1.52 | 0.131 | |
| Porosity × Absorption × Country (Ch.) | 0.48 | 0.22 | 2.21 | 0.028 | * |
| Porosity × Absorption × Country (Vt.) | 0.08 | 0.25 | 0.33 | 0.740 | |
| Porosity × Absorption × Religion (Ch.) | -0.03 | 0.09 | -0.27 | 0.785 | |
| Porosity × Country (Gh.) × Religion (Ch.) | 0.03 | 0.12 | 0.22 | 0.827 | |
| Porosity × Country (Th.) × Religion (Ch.) | -0.22 | 0.15 | -1.41 | 0.161 | |
| Porosity × Country (Ch.) × Religion (Ch.) | 0.28 | 0.16 | 1.81 | 0.071 | |
| Porosity × Country (Vt.) × Religion (Ch.) | -0.17 | 0.24 | -0.71 | 0.478 | |
| Absorption × Country (Gh.) × Religion (Ch.) | 0.14 | 0.11 | 1.21 | 0.227 | |
| Absorption × Country (Th.) × Religion (Ch.) | -0.09 | 0.14 | -0.65 | 0.518 | |
| Absorption × Country (Ch.) × Religion (Ch.) | -0.48 | 0.23 | -2.06 | 0.040 | * |
| Absorption × Country (Vt.) × Religion (Ch.) | 0.19 | 0.13 | 1.46 | 0.144 | |
| Porosity × Absorption × Country (Gh.) × Religion (Ch.) | 0.04 | 0.11 | 0.31 | 0.755 | |
| Porosity × Absorption × Country (Th.) × Religion (Ch.) | 0.25 | 0.18 | 1.37 | 0.171 | |
| Porosity × Absorption × Country (Ch.) × Religion (Ch.) | -0.62 | 0.22 | -2.88 | 0.004 | ** |
| Porosity × Absorption × Country (Vt.) × Religion (Ch.) | 0.31 | 0.25 | 1.22 | 0.224 |
r12 <- lm(spirit_score_std ~ por_score_std * abs_score_std * country * site, d1)
# summary(r12)
regtab_fun(r12) %>% regtab_style_fun(row_emph = 2:3)
r3 <- lmer(other_score_std ~ por_score_std * abs_score_std +
(1 | country/site/religion),
# (1 | country/religion),
data = d1)
summary(r3)
rsquared(r3)
regtab_fun(r3, por_name = "Porosity Vignettes") %>%
regtab_style_fun(row_emph = 2:3)
regtab_ran_fun(r3) %>% regtab_style_fun()
rr3 <- lmer(other_recoded_std ~ por_score_std * abs_score_std +
(1 | country/site/religion),
# (1 | country/religion),
data = d1r)
summary(rr3)
piecewiseSEM::rsquared(rr3)
regtab_fun(rr3, por_name = "Porosity Vignettes") %>%
regtab_style_fun(row_emph = 2:3)
regtab_ran_fun(rr3) %>% regtab_style_fun()